% msm_ul.m
% 
% Matlab code for
% 
% "Why has the U.S. economy stagnated since the Great Recession?"
% by Yunjong Eo and James Morley, Review of Economics and Statistics, forthcoming
% 
% Structural break dates (two breaks)
% 
% gbdate: trend break 
% bdate: volatility
% 
% Sample period: 1947Q2 - 2018Q4 after taking log difference


clear, clc;
warning('off');


%Real GDP 1947.1 - 2018.4
data = dlmread('rgdp2018q4.txt');
dates = 1947.25:0.25:2018.75; %report results from 1947Q2

yy = log(data(:,1));
y = 100*(yy(2:size(yy,1))-yy(1:size(yy,1)-1));
T = size(y,1);

%Globals
mstarU=5;              %mstarU = 7, 6, 5, 4, 3, 2, or 1 for U-shaped recession
mstarL=mstarU;         %mstarL = 7, 6, 5, 4, 3, 2, or 1 for partial BB w/ L-shaped recession; not used in benchmark, but allows for more general model
growth_bdate = 236;         % 2006Q1 break for trend growth 
vol_bdate = 149;            % 1984Q2 break for volatility      


%% MLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prmtr_in = -zeros(10, 1); % Initial parameter values

prmtr_in =[-3.5;-3.9;0.6;0.9;0.9;-1.3;-2.0;-0.4;-0.2;-1.7];


% Optimization routine.
fprintf('Optimization in progress...\n');
options=optimoptions('fminunc','MaxFunctionEvaluations',15000,'MaxIterations',10000,'FiniteDifferenceType','central');
[xout,fout,cout,out,gout,hout] = fminunc(@(prmtr)lik_fcn(prmtr,T,y,mstarU,mstarL,vol_bdate,growth_bdate),prmtr_in,options);
% Returns paramter estimates, -LL value, code, gradient, hessian

% Final parameter values
prm_fnl = trans(xout);

% Use Hessian to find parameter standard errors
hessn0 = -hout;
cov0_h = inv(-hessn0);
grdn_fnl = numeric_jacobian(@(xout)trans(xout),xout);
      
cov_h=grdn_fnl*cov0_h*grdn_fnl';
sd_fnl =sqrt(diag(cov_h)); %Standard errors of the estimated coefficients based on Hessian

results = [prm_fnl, sd_fnl];

% Creates output file to store results
output = fopen('msm_ul.txt','w');

% Print table of results into 'qmle_CMS_M.txt'
fprintf(output, 'log likelihood value is %1.8f\n', -fout);
fprintf(output, 'Parameter, QMLE parameter Estimates, Standard errors\n');
fprintf(output, 'p_01        %0.8f              %0.8f\n', results(1, :));
fprintf(output, 'p_02        %0.8f              %0.8f\n', results(2, :));
fprintf(output, 'p_11        %0.8f              %0.8f\n', results(3, :));
fprintf(output, 'p_22        %0.8f              %0.8f\n', results(4, :));
fprintf(output, 'mu0         %0.8f              %0.8f\n', results(5, :));
fprintf(output, 'mu1        %0.8f              %0.8f\n', results(6, :));
fprintf(output, 'mu2        %0.8f              %0.8f\n', results(7, :));
fprintf(output, 'mu00       %0.8f              %0.8f\n', results(8, :));
fprintf(output, 'sig_e1      %0.8f              %0.8f\n', results(9, :));
fprintf(output, 'sig_e2      %0.8f              %0.8f\n', results(10, :));
fclose(output);

  
prmtr = prm_fnl;

[pr_tt1,pr_tt2,pr_tt3,pr_tl1,pr_tl2,pr_tl3,yhatvec,gapvec] = filtered_inference(prmtr,T,y,mstarU,mstarL,vol_bdate,growth_bdate);



RDSS = yy(2:T+1) - gapvec/100;

[pr_sm1,pr_sm2,pr_sm3] = smoothed_probs(prmtr,T,pr_tt1,pr_tt2,pr_tt3,pr_tl1,pr_tl2,pr_tl3);   % get smoothed probabilities

pr_smC = pr_sm2+pr_sm3;

%=========================================================================%
% Figures
%=========================================================================%

figure

subplot(2,1,1); 
h1=NBERbc(dates,y,{'-'},3,{'r'});
hold on
h2=plot(dates,yhatvec,'-.','LineWidth',3);
hold on
plot([dates(1) dates(end)],zeros(2,1),'-k','LineWidth',2)
%ylim([-6 5])
xlabel("Quarter"); title("Output Growth");
legend([h1 h2],{'real GDP Growth (Quarterly)','Mean Growth'})
set(gca,'FontSize',16)
hold off

subplot(2,1,2);
h1=NBERbc(dates,gapvec,{'-'},3,{'r'});
hold on
plot([dates(1) dates(end)],zeros(2,1),'-k','LineWidth',2)
%ylim([-6 5])
xlabel("Quarter"); title("Output Gap");
set(gca,'FontSize',16)
hold off


figure
subplot(2,1,1);
h1=NBERbc(dates,pr_smC,{'-'},3,{'r'});
hold on
plot([dates(1) dates(end)],zeros(2,1),'-k','LineWidth',2)
%ylim([0 1])
xlabel("Quarter"); title("Probability of Contraction");
set(gca,'FontSize',16)
hold off

subplot(2,1,2); 
h1=NBERbc(dates,pr_sm2,{'-'},3,{'r'});
hold on
h2=plot(dates,pr_sm3,'-.','LineWidth',3);
hold on
plot([dates(1) dates(end)],zeros(2,1),'-k','LineWidth',2)
%ylim([0 1])
xlabel("Quarter"); title("Probabilities of Contraction Regimes");
legend([h1 h2],{'Prob of L','Prob of U'})
set(gca,'FontSize',16)
hold off





